Method and apparatus for evaluating polynomials and rational functions

ABSTRACT

Disclosed herein are a computer-processing method and apparatus for computing values of polynomials or rational functions. A mathematical software library can advantageously embody the concepts of this invention. The method can be adapted to compute values for non-elementary,special functions, for example ERF, ERFC, LGAMMA, and Bessel functions. The steps for polynomial evaluation include presenting input data that includes coefficients of polynomial p(x), x, a predetermined x i , and p(x i ), building polynomial c(x) having coefficients so that polynomial p(x) is expressible as: 
       p ( x )= p ( x   i )+{ x−x   i   }·c ( x ), 
     determining each coefficient of polynomial c(x), determining a value of polynomial c(x), and constructing a value of polynomial p(x) by determining: 
       p ( x )= p ( x   i )+{ x−x   i   }·c ( x ). 
     The method can be adapted for providing a value for a rational function 
       r ( x )= p ( x )/ q ( x ), 
     which is a ratio of a numerator polynomial p(x) and a denominator polynomial q(x).

FIELD OF THE INVENTION

[0001] The present invention relates to a method and apparatus forcomputing values of mathematical expressions, and more specifically to acomputer processing method and computer apparatus for computing binaryrepresentations of numerical values of polynomials and rationalfunctions.

BACKGROUND OF THE PRESENT INVENTION

[0002] Polynomials and rational functions, which are quotients ofpolynomials, are used, for example, by various branches of appliedscience for determining numerical values of mathematical expressions formodeling a property of a physical system, for example a rate of air flowover an airfoil surface. Polynomials can be used to approximatecomplicated mathematical expressions. Various methods, for exampleHomer's Rule that was disclosed in 1819, are used for computing valuesof polynomials, and provide a degree of accuracy that may not beacceptable in certain situations, depending on the particular polynomialand the accuracy required. G. E. Forsythe et al in “Computer Methods forMathematical Computations” Prentice-Hall (1977) discloses, in Section4.2, using Homer's Rule for computing a numerical value of a polynomial.

[0003] Steps for computing binary representations of numbers can createan unacceptably large deviation between a computed binary representationand its theoretical numerical value due to successive rounding errors.This can be an intolerable situation when a higher degree of accuracy isrequired. Various methods can improve computation accuracy but they mayrequire a significant increase in processing time and/or hardware.

[0004] Agarwal et al in “New Scalar and Vector Elementary Functions forthe IBM System/370” IBM Journal of Research and Development: Vol.30 No.2(March 1986), and Gal et al in “An Accurate Elementary MathematicalLibrary for the IEEE Floating Point Standard” ACM Transactions onMathematical Software Vol.17 No.1 (March 1991) pp. 26 to 45, appear todisclose methods that include a lookup table having non-equidistantlyspaced entries for computing values for elementary mathematicalfunctions.

[0005] Markstein in “Computation of Elementary Functions on the IBM RISCSystem/6000 Processor” IBM Journal of Research and Development Vol.34No.1 (January 1990) appears to disclose a method for computing values ofa function g(x) near a given point x_(i).

[0006] Gustavson et al in “Elementary Function Subroutines” Departmentof Mathematical Sciences of IBM Watson Research Center (January 1985)appears to disclose a program implementation for a software mathematicallibrary.

OBJECT OF THE PRESENT INVENTION

[0007] An object of the present invention is to obviate thedisadvantages of the prior art.

SUMMARY OF THE PRESENT INVENTION

[0008] A first aspect of the present invention provides amachine-processing method for computing a property of a mathematicallymodelled physical system, the steps including: reading, via a machineprocessing unit, input data including a value for each identifiedordered coefficient of a first polynomial p(x) representing theproperty, the polynomial p(x) being expressed as

p(x)=Σ(P _(j) ·x ^(j)) where j=0 to n,

[0009] a value of a quantity x, a value of a predetermined x_(i), and avalue of a predetermined p(x) correspondingly paired with saidpredetermined x_(i); building, via the machine processing unit, a valueof a second polynomial c(x) having ordered coefficients, the secondpolynomial c(x) being expressible as:

c(x)=Σ(C _(k) ·x ^(k)) where k=0 to (n−1)

[0010] so that the first polynomial p(x) is expressible as:

p(x)=p(x _(i))+{x−x _(i) }·c(x),

[0011] including the steps of: determining, via the machine processingunit, a value for each ordered coefficient of the second polynomial c(x)by generating a plurality of machine processing unit signals todetermine each of the ordered coefficients of the second polynomial c(x)from:

C _(k)=Σ(P _((k+1+j)) ·x _(j) ^(i)) where j=0 to (n−1−k);

[0012] determining, via the machine processing unit, a value of thesecond polynomial c(x) by generating a plurality of machine processingunit signals to determine:

c(x)=Σ(C _(k) ·x ^(k)) where k=0 to (n−1);

[0013] and constructing, via the machine processing unit, a value of thefirst polynomial p(x) by generating a plurality of machine processingunit signals to determine

p(x)=p(x _(i))+{x−x _(i) }·c(x);

[0014] and outputting, via the machine-processing unit, the value of thefirst polynomial p(x) representing said property of the mathematicallymodelled physical system.

[0015] Preferably the machine-implementable method is further adapted sothat a difference between x and x_(i) is sufficiently small to achieve adesired accuracy of a final computation of the machine representation ofa numerical value of the first polynomial p(x).

[0016] Preferably the machine-implementable method is further adapted sothat the step of reading the input data includes reading, via themachine processing unit, the input data from a machine-readable medium.

[0017] Preferably the machine-implementable method is further adapted sothat the ordered coefficients of the second polynomial c(x) are computedfrom a mathematical expression derivable from:

C _(k)=Σ(P _((k+1+j)) ·x _(i) ^(j)) where j=0 to (n−1−k).

[0018] Preferably the machine-implementable method is further adapted sothat the mathematical expression is a mathematical recurrenceexpression.

[0019] Preferably the machine-implementable method is further adapted sothat the mathematical recurrence expression is a forward mathematicalrecurrence expression.

[0020] Preferably the machine-implementable method is further adapted sothat the mathematical recurrence expression is a backward mathematicalrecurrence expression.

[0021] Preferably the machine-implementable method is further adapted toimplement the backward mathematical recurrence expression by includingfurther steps for: equating, via the machine-processing unit, a value ofa highest ordered coefficient of the second polynomial c(x) to a valueof an identified highest ordered coefficient of the first polynomialp(x) by generating a plurality of machine processing unit signals todetermine:

C _(n−1) =P _(n);

[0022] and computing, via a machine processing unit, a value for eachlower ordered coefficient of the second polynomial c(x) by generating aplurality of machine processing unit signals to determine:

C _(k−1) =x _(i) ·C _(k) +P _(k) where k=(n−1) to 1.

[0023] Preferably the machine-implementable method is further adapted sothat the predetermined x_(i) is selected from a set of predeterminedvalues of x_(i).

[0024] Preferably the machine-implementable method is further adapted sothat the predetermined x_(i) is a closest member of the set to theidentified x.

[0025] Preferably the machine-implementable method is further adapted sothat the step of determining a value of the second polynomial c(x) iscomputed by using Horner's Rule.

[0026] Preferably the machine-implementable method is further adaptedfor determining a value of a denominator polynomial q(x) havingidentified ordered denominator coefficients, the denominator polynomialq(x) being expressible as:

q(x)=Σ(i Q_(j) ·x ^(j)) where j=0 to m,

[0027] including further steps of: adapting the input data to furtherinclude a value for each identified ordered denominator coefficient ofthe denominator polynomial q(x), a value of a predetermined q(x)correspondingly paired with the predetermined x_(i); and determining,via the machine processing unit, a value of another polynomial d(x)having ordered denominator coefficients, the another polynomial d(x)being expressible as:

d(x)=Σ(D _(k) ·x ^(k)) where k=0 to (m−1)

[0028] so that the denominator polynomial q(x) is expressible as:

q(x)=q(x _(i))+{x−x _(i) }·d(x),

[0029] and a value for the denominator polynomial is resolved.

[0030] Preferably the machine-implementable method is further adapted sothat the first polynomialp(x) is a numerator polynomialp(x), andp(x)−p(x_(i)) is computed, and p(x_(i)) is not read.

[0031] Preferably the machine-implementable method is further adaptedfor determining a value of a rational function r(x) being expressible asa quotient of the numerator polynomial p(x) and the denominatorpolynomial q(x) expressed as ${{r(x)} = \frac{p(x)}{q(x)}},$

[0032] including further steps of: adapting the input data to furtherincluding a value of a predetermined r(x_(i)) correspondingly pairedwith the predetermined x_(i); constructing, via the machine processingunit, a value of the rational function r(x) by generating a plurality ofmachine processing unit signals to determine:${r(x)} = {{{r\left( x_{i} \right)} \cdot \left( {1 - \frac{{q(x)} - {q\left( x_{i} \right)}}{q(x)}} \right)} + {\frac{{p(x)} - {p\left( x_{i} \right)}}{q(x)}.}}$

[0033] Preferably the machine-implementable method is further adapted sothat the rational function is an approximation to a Bessel function.

[0034] Preferably the machine-implementable method is further adapted sothat the rational function is an approximation to an error function(ERF).

[0035] Preferably the machine-implementable method is further adapted sothat the rational function is an approximation to a complementary errorfunction (ERFC).

[0036] Preferably the machine-implementable method is further adapted sothat the rational function is an approximation to a log gamma function(LGAMMA).

[0037] Preferably the machine-implementable method is further adapted sothat all values are machine representations of some numerical value, themachine processing unit is a computer processing unit, and each machinerepresentation is a binary representation operable with the computerprocessing unit, and the machine-readable medium is a computer-readablemedium.

[0038] A second aspect of the present invention provides acomputer-readable program product having computer executableinstructions for instructing a computer, the instructions tangiblyembodying the machine-implementable method of the first aspect of thepresent invention.

[0039] A third aspect of the present invention provides acomputer-readable mathematical software routine library includingcomputer executable instructions for instructing a computer, theinstructions tangibly embodying the machine-implementable method of thefirst aspect of the present invention.

[0040] Preferably, the computer-readable mathematical software routinelibrary is further adapted so that the library is operatively associatedwith a software programming language.

[0041] A fourth aspect of the present invention provides a machine forcomputing a property of a mathematically modelled physical system, thesteps including: reading, via a machine processing unit, input dataincluding a value for each identified ordered coefficient of a firstpolynomial p(x) representing the property, the polynomial p(x) beingexpressed as

p(x)=Σ(P _(j) ·x ^(j)) where j=0 to n,

[0042] a value of a quantity x, a value of a predetermined x_(i), and avalue of a predetermined p(x_(i)) correspondingly paired with thepredetermined x_(i); building, via the machine processing unit, a valueof a second polynomial c(x) having ordered coefficients, the secondpolynomial c(x) being expressible as:

c(x)=Σ(C _(k) ·x ^(k))where k=0 to (n−1)

[0043] so that the first polynomial p(x) is expressible as:

p(x)=p(x _(i))+{x−x _(i) }·c(x),

[0044] including the steps of: determining, via the machine processingunit, a value for each ordered coefficient of the second polynomial c(x)by generating a plurality of machine processing unit signals todetermine each of the ordered coefficients of the second polynomial c(x)from:

C _(k)=Σ(P _((k+1+j)) ·x _(j) ^(i)) where j=0 to (n−1−k);

[0045] determining, via the machine processing unit, a value of thesecond polynomial c(x) by generating a plurality of machine processingunit signals to determine:

c(x)=Σ(C _(k)·x^(k)) where k=0 to (n−1);

[0046] constructing, via the machine processing unit, a value of thefirst polynomial p(x) by generating a plurality of machine processingunit signals to determine :

p(x)=p(x _(i))+{x−x _(i) }·c(x);

[0047] and outputting, via the machine-processing unit, the value of thefirst polynomial p(x) representing the property of the mathematicallymodelled physical system.

[0048] Preferably, the machine is further adapted so that a differencebetween x and x_(i) is sufficiently small to achieve a desired accuracyof a final computation of the machine representation of a numericalvalue of the first polynomial p(x).

[0049] Preferably, the machine is further adapted so that the means forreading the input data includes means for reading, via the machineprocessing unit, the input data from a machine-readable medium.

[0050] Preferably, the machine is further adapted so that the orderedcoefficients of the second polynomial c(x) are computed from amathematical expression derivable from:

C _(k)=Σ(P _((k+1+j)) ·x _(i) ^(j)) where j=0 to (n−1−k).

[0051] Preferably, the machine is further adapted so that themathematical expression is a mathematical recurrence expression.

[0052] Preferably, the machine is further adapted so that themathematical recurrence expression is a forward mathematical recurrenceexpression.

[0053] Preferably, the machine is further adapted so that themathematical recurrence expression is a backward mathematical recurrenceexpression.

[0054] Preferably, the machine is further adapted to implement thebackward mathematical recurrence expression by further including: meansfor equating, via the machine processing unit, a value of a highestordered coefficient of the second polynomial c(x) to a value of anidentified highest ordered coefficient of the first polynomial p(x) bygenerating a plurality of machine processing unit signals to determine:

C _(n−1) =P _(n);

[0055] and means for computing, via the machine processing unit, a valuefor each lower ordered coefficient of the second polynomial c(x) bygenerating a plurality of machine processing unit signals to determine:

C _(k−1) =x _(i) ·C _(k) +P _(k) where k=(n−1) to 1.

[0056] Preferably, the machine is further adapted so that thepredetermined x_(i) is selected from a set of predetermined values ofx_(i).

[0057] Preferably, the machine is further adapted so that thepredetermined x_(i) is a closest member of the set to the identified x.

[0058] Preferably, the machine is further adapted for determining meansfor determining a value of the second polynomial c(x) is computed byusing Homer's Rule.

[0059] Preferably, the machine is further adapted for determining avalue of a denominator polynomial q(x) having identified ordereddenominator coefficients, the denominator polynomial q(x) beingexpressible as:

q(x)=Σ(Q _(j) ·x ^(j)) where j=0 to m,

[0060] including further steps of: adapting the input data to furtherinclude a value for each identified ordered denominator coefficient ofthe denominator polynomial q(x), and a value of a predetermined q(x_(i))correspondingly paired with the predetermined x_(i); and determining,via the machine processing unit, a value of another polynomial d(x)having ordered denominator coefficients, the another polynomial d(x)being expressible as:

d(x)=Σ(D _(k) ·x ^(k)) where k=0 to (m−1)

[0061] so that the denominator polynomial q(x) is expressible as:

q(x)=q(x _(i))+{x−x _(i) }·d(x),

[0062] and a value for the denominator polynomial is resolved.

[0063] Preferably, the machine is further adapted so that the firstpolynomial p(x) is a numerator polynomial p(x), and p(x)−p(x_(i)) iscomputed, and p(x_(i)) is not read.

[0064] Preferably, the machine is further adapted for determining avalue of a rational function r(x) being expressible as a quotient of thenumerator polynomial p(x) and the denominator polynomial q(x) expressedas ${{r(x)} = \frac{p(x)}{q(x)}},$

[0065] including further steps of: adapting the input data to furtherincluding a value of a predetermined r(x_(i)) correspondingly pairedwith the predetermined x_(i); and constructing, via the machineprocessing unit, a value of the rational function r(x) by generating aplurality of machine processing unit signals to determine:${r(x)} = {{{r\left( x_{i} \right)} \cdot \left( {1 - \frac{{q(x)} - {q\left( x_{i} \right)}}{q(x)}} \right)} + {\frac{{p(x)} - {p\left( x_{i} \right)}}{q(x)}.}}$

[0066] Preferably, the machine is further adapted so that the rationalfunction is an approximation to a Bessel function.

[0067] Preferably, the machine is further adapted so that the rationalfunction is an approximation to an error function (ERF).

[0068] Preferably, the machine is further adapted so that the rationalfunction is an approximation to a complementary error function (ERFC).

[0069] Preferably, the machine is further adapted so that the rationalfunction is an approximation to a log gamma function (LGAMMA).

[0070] Preferably, the machine is further adapted so that all values aremachine representations of some numerical value, the machine processingunit is a computer processing unit, each machine representation is abinary representation operable with the computer processing unit, andthe machine-readable medium is a computer-readable medium.

[0071] A fifth aspect of the present invention provides a machine havinga computer-readable program product having computer executableinstructions for instructing a computer to embody the machine of thefourth aspect of the present invention.

[0072] A sixth aspect of the present invention provides a machine havinga computer-readable mathematical software routine library includingcomputer executable instructions for instructing a computer to embodythe machine of the fourth aspect of the present invention.

[0073] Preferably the computer-readable mathematical software routinelibrary of the sixth aspect of the present invention is further adaptedso that the library is operatively associated with a softwareprogramming language.

DETAILED DESCRIPTION OF THE FIGURES OF THE PRESENT INVENTION

[0074] To illustrate aspects of the present invention, the followingfigures are used, in which:

[0075]FIG. 1 depicts a programming flow chart for computing a value of apolynomial p(x) in accordance with the present invention;

[0076]FIG. 2 depicts a programming flow chart for computing a value of arational function p(x)/q(x) in accordance with the present invention;and

[0077]FIG. 3 depicts a programming flow chart for using a subroutine forcomputing a value of a Bessel function J₀(x) in accordance with thepresent invention.

DETAILED DESCRIPTION OF A PREFERRED EMBODIMENT OF THE INVENTION

[0078] The present invention will be described with reference to anexemplary context of a computer processing method and apparatus forcomputing a binary representation of a numerical value of a polynomialp(x).

[0079] For embodiments of the present invention, references to computingor calculating values or numbers will be understood to refer toprogramming instructions for instructing a computer to compute binaryrepresentations of numerical values of a mathematical expression orvariable. It is understood that computing machinery can manipulate andtransform binary representations of numerical values or numbers.

[0080] A polynomial p(x) can be mathematically expressed as shown in thefollowing equation:

p(x)=Σ P _(j) x ^(j) where j=0 to n in which p(x)=P ₀ +P ₁ ·x+P ₂ ·x ² +. . . +P _(n) ·x ^(n)  Equ. 1

[0081] Ordered coefficients of polynomial p(x) are P₀, . . . , P_(n).Polynomial p(x) has n+1 ordered coefficients and is said to have degreen. A numerical value of polynomial p(x) can be computed by reading eachcoefficient of polynomial p(x) and an identified x and using steps inEquation 1. Following Equation 1 directly would require a computerprocessor to perform $\frac{n\left( {n + 1} \right)}{2}$

[0082] multiplications and n additions. To significantly reduce theprocessing time, Homer's Rule can be applied to Equation 1, in whichpolynomial p(x) can be expressed in Equation 1a, as shown below:

p(x)=P ₀ +x(P ₁ +x(P ₂ +x( . . . (P _(n−1) +x·P _(n)) . . . )))  Equ. 1a

[0083] By following Equation 1a, a program can perform n multiplicationsand n additions to compute a number for polynomial p(x) using lessprocessing time than a program that follows Equation 1.

[0084] A polynomial p(x) at x=x_(i) can be expressed as follows:

p(x _(i))=ΣP _(j) x _(i) ^(j) where j=0 to n  Equ. 2

[0085] Therefore, a difference between Equation 1 and Equation 2 can beexpressed as follows:

p(x)−p(x _(i))=ΣP _(j)(x ^(j) −x _(i) ^(j)) where j=0 to n  Equ. 3

[0086] Expanding Equation 3 can provide the following expression:

p(x)−p(x _(i))=Σ{P _(j)(x−x _(i))·Σ[x ^(k) ·x _(i) ^((j−1−k))]} wherej=0 to n, k=0 to (j−1)  Equ. 4

[0087] For Equation 4, in the first sum j=0 to n, and in the second sumk=0 to (j−1). It is appreciated that Equation 4 can be expressed asfollows:

p(x)−p(x _(i))=(x−x _(i))·Σ{(ΣP _((k+1+j)) ·x _(i) ^(j))·x ^(k)} wherek=0 to (n−1), j=0 to (n−1−k)  Equ. 4a

[0088] For Equation 4a, in the first sum k=0 to (n−1), and in the secondsum j=0 to (n−1−k) By following Equation 4a, a value for polynomial p(x)can be computed by reading coefficients of the polynomial (P₁, P₂, . . ., P_(n)), and x, x_(i), and p(x_(i)). Equation 4a can be expressed as acombination of Equation 5 and Equation 6,as follows:

p(x)−p(x _(i))={x−x _(i) ·}c(x)={x−x _(i)}·Σ(C _(k) ·x ^(k)) where k=0to (n−1)  Equ. 5

C _(k) =ΣP _((k+1+j)) ·x _(i) ^(j) where j=0 to (n−1−k)  Equ. 6

[0089] Equation 5 includes an expression for a new polynomial c(x)having coefficients, from C₀, . . . , C_((n−1)), that are not known apriori and therefore need to be determined. Equation 6 can be used fordetermining values of the coefficients of new polynomial c(x) byinvolving the coefficients of polynomial p(x) and an identified x_(i).Optionally, other expressions for determining values for eachcoefficient of new polynomial c(x) can be derived from Equation 6.

[0090] It is appreciated that Homer's Rule could be used to compute avalue for polynomial c(x) after coefficients of polynomial c(x) havebeen computed by using Equation 6 or optionally from another expressionthat could be derived from Equation 6. By using Homer's Rule, Equation 5can be expressed as follows:

p(x)−p(x _(i))={x−x _(i) }·{C ₀ +x(C ₁ +x(C ₂ +x( . . . (C _(n−2) +x·C_(n−1)) . . . )))}  Equ. 5a

[0091] Processing steps that follow Equation 5a begin within theinnermost set of parentheses and proceed outwards to the outermostbrackets. Calculated coefficients of new polynomial c(x) cansubsequently be used in Equation 5a for computing a numerical value ofpolynomialp(x).

[0092] Optionally, unknown coefficients of new polynomial c(x) can becomputed by applying the principle of mathematical recurrence onEquation 6. For example, a forward mathematical recurrence can allowcomputation of coefficients of polynomial c(x) by beginning with thelowest ordered coefficient, C₀, and proceeding to the highest orderedcoefficient, C_((n−1)). Optionally, a backward mathematical recurrencecan compute coefficients of polynomial c(x) by beginning with thehighest ordered coefficient, C_((n−1)) and proceeding to the lowestordered coefficient, C₀. With each step of backwardly computing acoefficient of polynomial c(x), a next outer most bracketed term ofEquation 5a can be advantageously computed.

[0093] Optionally, a backward recurrence for determining eachcoefficient of new polynomial c(x) can be realized by using Equation 7and Equation 8. Each subsequent lower-ordered unknown coefficient of newpolynomial c(x) can be computed by using Equation 8.

C _(n−1) =P _(n)  Equ. 7

C _(k−1) =x _(i) ·C _(k) +P _(k) where k=(n−1) to 1  Equ. 8

[0094] To compute a numerical value of a polynomial p(x), a computer canread input data that includes a value for each identified orderedcoefficient of the first polynomial p(x), a value of a quantity x, avalue of a predetermined x_(i), a value of a predetermined p(x_(i)) thatis correspondingly paired with said predetermined x_(i). An optionalstep can select a nearest x_(i) that is closest to an identified x, inwhich the nearest x_(i) is selected from a set of predetermined x_(i),and a nearest predetermined p(x_(i)) corresponds to the nearest x_(i).Optionally, coefficients of polynomial c(x) can be computed by followingEquation 7 and Equation 8. A value for polynomial p(x) can then becomputed by following Equation 5a. The operations of computing thecoefficients of polynomial c(x) in Equation 8 can be interleaved withcomputing the nested parenthesis in Equation 5a so that each time a newC_(k) is computed by Equation 8 the next outer nested parenthesis ofEquation 5a is computed next.

[0095] Provided p(x)−p(x_(i)) does not have a multiple root at x=x_(i),the relative round-off error in computing p(x)−p(x_(i)) with Equation 5can be kept small by keeping x sufficiently close to x_(i). This can beachieved by having available a number of sufficiently closely spacedx_(i)'s covering the range of x to be handled, and using the nearestx_(i) to x. To obtain a reasonable result that would be satisfactory fora given computational purpose, a difference between x and x_(i) can besufficiently small enough to achieve a desired accuracy of a computedvalue of polynomial p(x).

[0096] Optionally, accuracy can be further improved by choosing x_(i)such that p(x_(i)) has extra accuracy beyond the precision of thefloating-point number system of the computer. This can be achieved byselecting x_(i) to be a floating-point number, by exhaustively searchingfrom a starting desired value for x_(i), until an x_(i) is found forwhich p(x_(i)) is very close to a floating-point number (i.e. closerthan half the distance between adjacent floating-point numbers). Such afloating-point approximation to p(x_(i)) has an accuracy beyond theprecision of the floating-point number system.

[0097]FIG. 1 shows the preferred embodiment of the present invention,which is represented as a computer processing method embodied in aprogramming flow chart for computing a machine representation of anumerical value of a polynomial p(x). Optionally, coefficients of p(x)can be placed in a vector of length n. The method optionally includes astep of accessing a table having a set of predetermined values of x_(i)and a corresponding set of predetermined values of p(x_(i)) each pairedwith a corresponding x_(i). A table can contain values for the followingexpressions:

x _(i) , p(x _(i)) for i=1 to N.  Equ. 8a

[0098] In step 10, the method of the preferred embodiment begins. Instep 15, each coefficient of polynomial p(x) and an identified x can beread by a computer processing unit. In step 20, the computer optionallyselects a predetermined x_(i) from a set of predetermined values ofx_(i). Preferably, the predetermined x_(i)is a closest member of the setto the identified x. In step 30, the initial conditions for aprogramming loop are set. Variable ksum is an ongoing summation thatrepresents the sum in Equation 5a. Variable ksum is initially set equalto coefficient P_(n) of polynomial p(x). Variable ck is a temporaryvalue for any coefficients of polynomial c(x), in which variable ck iscomputed to be a specific coefficient of polynomial c(x) that isdetermined in step 50 as will be shown below. Variable ck is initiallyset equal to P_(n), which is the highest ordered coefficient ofpolynomial p(x). Variable k is an index counter for a processing loop,in which variable k is initially set equal to (n−1), which is one lessthan the degree of polynomial p(x). In step 40, a query is performed todetermine whether a value for variable k is equal to zero. If thenumerical value for variable k is greater than zero, the method proceedsto step 50. If the numerical value for variable k is equal to zero, themethod proceeds to step 60. In step 50, a value of variable ck iscomputed, which is coefficient C_(k) by using the backward recurrence ofEquation 8 that involves a higher ordered coefficient of polynomialp(x), which is P_(k), and a previously determined coefficient, C_(k−1),of polynomial c(x). Next, an updated value for variable ksum iscomputed, which is an ongoing summation of the sum in Equation 5a, bydetermining a sum for the next outer bracketed term of Equation 5a. Forremaining iterations of step 50, the summation represented in Equation5a progresses from an innermost set of parentheses and proceedsoutwardly for each iterative step in an ongoing manner towards theoutermost brackets. An updated value of variable k is computed so thatthe sum in Equation 5a can be resolved towards its outermost brackets.In step 60, the method uses a value for polynomial c(x), which isexpressed in FIG. 1 as variable ksum after variable k equals zero instep 40. A value of polynomial p(x) is computed as a rearrangement ofEquation 5a, which is expressed as

p(x)=p(x _(i))+(x−x _(i))*c(x).

[0099] In step 70, the method for the preferred embodiment can stop orcan be repeated from step 10.

[0100] A second embodiment of the present invention provides a computerprocessing method for computing a machine representation of a numericalvalue of a rational function r(x). The rational function r(x) isexpressed as a quotient of a numerator polynomial p(x) having a set of npredetermined coefficients from P₀, P₁, . . . , P_(n), so that

p(x)=ΣP _(j) ·x ^(j) where j=0 to n,

[0101] and a denominator polynomial q(x) having m predeterminedcoefficients from Q₀, Q₁, . . . , Q_(m), so that

q(x)=ΣQ _(j) ·x ^(j) where j=0 to m.

[0102] A rational function r(x) can be defined in the followingexpression: $\begin{matrix}{{r(x)} = \frac{p(x)}{q(x)}} & \text{Equ.~~8b}\end{matrix}$

[0103] From Equation 8b, it follows that: $\begin{matrix}{{r(x)} = {{{r\left( x_{i} \right)} \cdot \left( {1 - \frac{{q(x)} - {q\left( x_{i} \right)}}{q(x)}} \right)} + \frac{{p(x)} - {p\left( x_{i} \right)}}{q(x)}}} & {{Equ}.\quad 9}\end{matrix}$

[0104] The second method can optionally include access to a table havinga set of predetermined values of x_(i) and a corresponding set ofpredetermined values of q(x_(i)) each paired with a corresponding x_(i),and a corresponding set of predetermined values of r(x_(i)) each pairedwith a corresponding x_(i). It is understood that each unique x_(i) ispaired with a correspondingly specific q(x_(i), and a correspondinglyspecific r(x_(i)). A spacing for table values of x_(i) can be set suchthat a resulting absolute value of (x−x_(i)) is small enough to obtain adesired accuracy. The x_(i) spacing required for a given desiredaccuracy is typically estimated with a perturbation analysis (usingTaylor series, for example), and verified by numerical testing. Thecloseness of x to x_(i) needed to achieve a given accuracy in thecomputed r(x) depends on the particular rational function r(x).

[0105] A table can contain values for the following expressions:

x _(i) , r(x _(i)),q(x _(i)) for i=1 to N  Equ. 10

r(x _(i)) and q(x _(i)) for i=1 to N  Equ. 11

[0106] Accuracy can be further improved by choosing the x_(i)'s suchthat the following terms shown in Equation 11 can have extra accuracybeyond the precision of the floating-point number system of a computer.This can be achieved by selecting x_(i) to be a floating-point number,by exhaustively searching from a starting desired value for x_(i), untilan x_(i) is found for which so that both the quantities in Equation 11are very close to floating-point numbers (closer than half the distancebetween adjacent floating-point numbers). Such floating-pointapproximations to r(x_(i)) and q(x_(i)) have an accuracy beyond theprecision of the floating-point number system.

[0107] The second method computes values of a numerator polynomial p(x)and a denominator polynomial q(x) by following Equation 5, as follows:

p(x)−p(x_(i))={x−x _(i) }·c(x)  Equ. 12a

q(x)−q(x _(i))={x−x _(i) }·d(x)  Equ. 12

[0108] Polynomial c(x) has (n−1) unknown coefficients and polynomiald(x) has (m−1) unknown coefficients. The second method can determinecoefficients of c(x) and d(x) and then determine a value for rationalfunction r(x).

[0109]FIG. 2 shows the second embodiment of the present invention, whichis represented as a computer processing method embodied in a programmingflow chart for computing a computer representation of a numerical valueof a rational function r(x)=p(x)/q(x). Optionally, steps for accessing atable having a set of various predetermined values for x_(i) a set ofvarious predetermined values for q(x_(i)), and a set of variouspredetermined values for r(x_(i)) can be included.

[0110] In step 110, the second method begins. In step 115, a computerreads values of each coefficient of polynomial p(x) and each coefficientof polynomial q(x), and an identified x. In step 120, a predeterminedx_(i) from a set of predetermined values of x_(i) can optionally beselected. Preferably, the predetermined x_(i) is a closest member of theset to the identified x. These predetermined values can be optionallystored in a table or they can be supplied to the computer in some othermeans. Optionally, the table includes values for various xi, variousq(xi) associated with the corresponding x_(i), and various r(x_(i))associated with the corresponding x_(i). Optionally, specific values forthe x_(i)'s, q(x_(i))'s, and r(x_(i))'s are chosen so that theq(x_(i))'s and r(x_(i))'s have extra accuracy beyond machine precisionin the manner of the accurate table method that is disclosed by theprior art. It is appreciated that the computed result will be moreaccurate if the table has extra accuracy. Step 130, step 140, and step150 are similar to step 30, step 40, and step 50 of FIG. 1, in that step130, step 140, and step 150 provide a value for variable ksum to be usedin step 160 as will be shown next. In step 160, a value for variable dp(i.e., delta p) is computed,which is a computed value for p(x)−p(x_(i))for a numerator polynomialp(x). Variable dp will be used in step 205 aswill be shown later. Step 170, step 180, and step 190 are similar tostep 130, step 140, and step 150 of FIG. 2, in that step 170, step 180,and step 190 provide a value for variable ksum to be used in step 200 aswill be shown next. In step 200 a value for variable dq (i.e., delta q),which is a computed value for q(x)−(x_(i)) for a denominator polynomialq(x) is computed. Variable dq will be used in step 202 and 205 as shownlater. In step 202 a numerical value for q(x) is computed by rearrangingEquation 12. In step 205 a value for a rational function r(x) iscomputed by following Equation 9. In step 210 the second method can stopor can be restarted from step 110.

[0111] A third embodiment of the present invention provides a computerprocessing method for computing binary representations of a value for aBessel function. A programming flow chart can be developed by adaptingFIG. 2 for computing a value for a Bessel function J₀(x), which is 20 aBessel function of the first kind of order 0. The third method can befurther adapted to compute a value for various Bessel functions, forexample J₀, J₁, Y₀, and Y₁, over a range of x (e.g. [0,8]) by firstapproximating the Bessel function by a piece-wise rational function.Bessel functions can typically be approximated by rational functionshaving numerator and denominator polynomials with the same degree (i.e.same number of coefficients).

[0112] To compute a value of a Bessel function J₀(x) using the thirdmethod, the first task would be to determine a range that an identifiedx belongs in. The following is a list of steps for determining theproper range to which x belongs: if the identified x is a NAN (i.e., nota legitimate number, such as $\frac{0}{0}$

[0113] or +/− INF (i.e., any non-zero number divided by zero), themethod does not compute a value for the Bessel function; if theidentified x is less than zero, the method makes the sign of theidentified x positive; if the identified x is greater than or equal tozero and less than 3.72*10⁻⁹, then the value of the Bessel function is1, assuming the result is returned in IEEE double-precision floatingpoint; if the identified x is greater than or equal to 3.72*10⁻⁹ andless than 1, then the method uses a Taylor series to compute a value forthe Bessel function; if the identified x is greater than or equal to 1and less than 1.5, then the method calls a subroutine that incorporatesthe present invention for determining a value for the Bessel function.The subroutine will be explained in FIG. 3. In each case where thesubroutine is called, arguments P, Q, T, and n of the subroutine areread along with values of x and x_(i) associated with the particularrange that x belongs with; if the identified x is greater than or equalto 1.5 and less than 1.9, then the method calls the subroutine fordetermining a value for the Bessel function; if the identified x isgreater than or equal to 1.9 and less than 2, then the method calls thesubroutine for determining a value for the Bessel function; if theidentified x is greater than or equal to 2 and less than 4, then themethod performs steps that are beyond the scope of the presentinvention. The steps involve accurate evaluation of a rationalapproximation to J₀(x)/(x−z₀), where z₀ is the first zero of J₀; if theidentified x is greater than or equal to 4 and less than 4.5, then themethod calls the subroutine for determining a value for the Besselfunction; if the identified x is greater than or equal to 4.5 and lessthan 5, then the method calls the subroutine for determining a value forthe Bessel function; if the identified x is greater than or equal to 5and less than 6, then the method performs steps that are beyond thescope of the present invention. The steps involve accurate evaluation ofa rational approximation to J₀(x)/(x−z₁), where z₁ is the first zero ofJ₀; if the identified x is greater than or equal to 6 and 6.4, then themethod calls the subroutine for determining a value for the Besselfunction; if the identified x is grater than 6.4 and less than 7, thenthe method calls the subroutine for determining a value for the Besselfunction; if the identified x is greater than or equal to 7 and lessthan 7.5, then the method calls the subroutine for determining a valuefor the Bessel function; if the identified x is greater than or equal to7.5 and less than 8, then the method calls the subroutine fordetermining a value for the Bessel function; and if the identified x isgreater than or equal to 8 then the method uses an asymptoticapproximation for determining a value for the Bessel function, which isbeyond the scope of the present invention.

[0114]FIG. 3 shows the third embodiment of the present invention, whichis represented as a computer processing method embodied in a programmingflow chart for computing a machine representation of a numerical valueof a Bessel function J₀(x).

[0115] In step 360 a subroutine, which can be called eval_rat, of theflow chart of the third embodiment begins. The subroutine inputs anargument x, at which the value of a rational function is desired to becomputed, and P and Q, which are vectors containing the coefficients ofthe numerator and denominator polynomials p(x) and q(x), respectively,of the rational function, and T, which is a table containing a set ofpredetermined x_(i)'s and corresponding r(x_(i))'s and q(x_(i))'s, andn, which is the degree of p(x) and q(x). In step 362 a computerprocessing unit reads values for each coefficient of numeratorpolynomial p(x) and denominator polynomial q(x), in which bothpolynomials have the same number or coefficients, and also reads anidentified x. In step 364, to achieve improved performance, a nearestx_(i) and other required data are optionally read from a table. In step366 various variables are initialized for subsequent use in aninstruction loop. One software loop is used because the numerator andthe denominator polynomials have the same degree (that is, share thesame number of coefficients). This loop is similar to a combination ofthe two loops involving steps 130 to 150 and 170 to 190 in FIG. 2. Instep 368, a test determines whether the loop has been completed. In step370, variables are computed in a manner similar to that combining steps150 and 190 in FIG. 1. In step 372, variables are computed in a mannersimilar to that combining steps 160 and 200 to 205 shown in FIG. 1. Step374 marks the end of the subroutine of the flow chart of the thirdembodiment.

[0116] For non-elementary special functions, for example the errorfunction (ERF), the complementary error function (ERFC), or the Besselfunctions,range-reduction properties, such as those available for theelementary functions, are not known. To directly apply prior arttable-driven techniques, for example those used for the elementaryfunctions, for these non-elementary special functions, would require acomputer processing unit to read coefficients of each of a large numberof polynomials, one associated with each table point. This can make thetable unacceptably large.

[0117] Computing steps that embody the methods of the present inventioncan be used to determine values for polynomials or rational functionsapproximating the non-elementary special function. Optionally,coefficients of a polynomial or rational function associated with eachindividual table point, x_(i), do not have to be stored; instead, thesecoefficients can be computed using simple formulas from only one (forpolynomials) or two (for rational functions) stored function values pertable point, resulting in tables of manageable size. In this sense, thepresent invention provides a benefit similar to that of range reductionformulas, for functions for which range reduction formulas do not exist.

[0118] Mathematical software libraries, including those associated witha software programming language or operating system, can be adapted toinstruct a general purpose digital computer to embody the presentinvention. It is appreciated that a software does not necessarily haveto belong to a programming language or operating system. There are alsostand-alone libraries of mathematical routines, for example, the IBM™Engineering and Scientific Subroutine Library, ESSL™.

[0119] There are many examples that illustrate the practical applicationof the present invention for mathematically modelling a property of aphysical system. C. M. Close in “The Analysis of Linear Circuits”,Harcourt, Brace & World (1966), discloses in Chapters 6 and 11 the useof polynomials for determining frequency response characteristics ofelectrical systems. The same concepts can be applied to non-electricalsystems, such as mechanical, hydraulic, acoustical, and thermal systems.R. C. Weyrick in “Fundamentals of Automatic Control”, McGraw-Hill BookCompany (1975), discloses in Chapters 2 and 6 the use of polynomials formathematically modelling physical systems for predicting a behavior of aphysical system. D. Halliday and R. Resnick in “Physics: Parts 1 and 2”,John Wiley & Sons (1978) disclose in Chapters 3 and 4 the use ofpolynomials for mathematically modelling and predicting spatialco-ordinates of impact of a projectile being released from an aeroplane.

[0120] The present invention can be implemented in a computer-readablecomputer memory program having a medium for storing and carrying theprogram to a general purpose computer for instructing a computerprocessing unit to implement the embodiments of the present invention. Acomputer-memory product can be a floppy disk or compact disk that can beadapted to carry software that tangibly embodies the present invention.

[0121] A computer processing method that tangibly embodies the presentinvention can be used for computing a number for a mathematical modeland for comparing the computed number against an observed behaviour of aphysical phenomenon.

[0122] The concepts of the present invention can be further extended toa variety of other applications that are clearly within the scope ofthis invention.

[0123] Having thus described the present invention with respect to apreferred embodiment as implemented, it will be apparent to thoseskilled in the art that many modifications and enhancements are possibleto the present invention without departing from the basic concepts asdescribed in the preferred embodiment of the present invention.Therefore, what is intended to be protected by way of letters patentshould be limited only by the scope of the following claims.

The embodiments of the invention in which an exclusive property orprivilege is claimed are defined as follows:
 1. A machine-processingmethod for computing a property of a mathematically modelled physicalsystem, the steps comprising: a) reading, via a machine processing unit,input data including a value for each identified ordered coefficient ofa first polynomial p(x) representing said property, said polynomial p(x)being expressed as p(x)=Σ(P _(j) ·x ^(j)) where j=0 to n,  a value of aquantity x, a value of a predetermined x_(i), and a value of apredetermined p(x_(i)) correspondingly paired with said predeterminedx_(i); b) building, via said machine processing unit, a value of asecond polynomial c(x) having ordered coefficients, said secondpolynomial c(x) being expressible as: C(x)=Σ(C _(k) ·x ^(k)) where k=0to (n−1)  so that said first polynomial p(x) is expressible as: p(x)=p(x_(i))+{x−x _(i) }·c(x),  comprising the steps of: i) determining, viasaid machine processing unit, a value for each ordered coefficient ofsaid second polynomial c(x) by generating a plurality of machineprocessing unit signals to determine each said ordered coefficient ofsaid second polynomial c(x) from: C _(k)=Σ(P _((k+1+j)) ·x _(i) ^(j))where j=0 to (n−1−k); ii) determining, via said machine processing unit,a value of said second polynomial c(x) by generating a plurality ofmachine processing unit signals to determine: C(x)=Σ(C _(k) ·x ^(k))where k=0 to (n−1); c) constructing, via said machine processing unit, avalue of said first polynomial p(x) by generating a plurality of machineprocessing unit signals to determine: p(x)=p(x _(i))+{x−x _(i) }·c(x); and d) outputting, via said machine-processing unit, said value of thefirst polynomial p(x) representing said property of the mathematicallymodelled physical system.
 2. The machine-implementable method of claim1, wherein a difference between x and x_(i) is sufficiently small toachieve a desired accuracy of a final computation of said machinerepresentation of a numerical value of said first polynomial p(x). 3.The machine-implementable method of claim 2 wherein the step of readingsaid input data comprises reading, via said machine processing unit,said input data from a machine-readable medium.
 4. Themachine-implementable method of claim 3 wherein said orderedcoefficients of said second polynomial c(x) are computed from amathematical expression derivable from: C _(k)=Σ(P _((k+1+j)) ·x _(i)^(j)) where j=0 to (n−1−k).
 5. The machine-implementable method of claim4 wherein said mathematical expression is a mathematical recurrenceexpression.
 6. The machine-implementable method of claim 5 wherein saidmathematical recurrence expression is a forward mathematical recurrenceexpression.
 7. The machine-implementable method of claim 5 wherein saidmathematical recurrence expression is a backward mathematical recurrenceexpression.
 8. The machine-implementable method of claim 7 furtheradapted to implement said backward mathematical recurrence expression bycomprising further steps for: e) equating, via said machine-processingunit, a value of a highest ordered coefficient of said second polynomialc(x) to a value of an identified highest ordered coefficient of saidfirst polynomial p(x) by generating a plurality of machine processingunit signals to determine: C _(n−1) =P _(n);  and f) computing, via amachine processing unit, a value for each lower ordered coefficient ofsaid second polynomial c(x) by generating a plurality of machineprocessing unit signals to determine: C _(k−1) =x _(i) ·C _(k) +P _(k)where k=(n−1) to
 1. 9. The machine-implementable method of claim 8wherein said predetermined x_(i) is selected from a set of predeterminedvalues of x_(i).
 10. The machine-implementable method of claim 9 whereinsaid predetermined x_(i) is a closest member of said set to saididentified x.
 11. The machine-implementable method of claim 10 whereinsaid step of determining a value of said second polynomial c(x) iscomputed by using Homer's Rule.
 12. The machine-implementable method ofclaim 11 for determining a value of a denominator polynomial q(x) havingidentified ordered denominator coefficients, said denominator polynomialq(x) being expressible as: q(x)=Σ(Q _(j) ·x ^(j)) where j=0 to m,comprising further steps of: g) adapting said input data to furtherinclude a value for each identified ordered denominator coefficient ofsaid denominator polynomial q(x), a value of a predetermined q(x_(i))correspondingly paired with said predetermined x_(i); and h)determining, via said machine processing unit, a value of anotherpolynomial d(x) having ordered denominator coefficients, said anotherpolynomial d(x) being expressible as: d(x)=Σ(D _(k) ·x ^(k)) where k=0to (m−1)  so that said denominator polynomial q(x) is expressible as:q(x)=q(x_(i))+{x−x _(i) }·d(x),  and a value for the said denominatorpolynomial is resolved.
 13. The machine-implementable method of claim 12wherein the first polynomial p(x) is a numerator polynomial p(x), andp(x)−p(x_(i)) is computed, and p(x_(i)) is not read.
 14. Themachine-implementable method of claim 13 for determining a value of arational function r(x) being expressible as a quotient of said numeratorpolynomial p(x) and said denominator polynomial q(x) expressed as${{r(x)} = \frac{p(x)}{q(x)}},$

comprising further steps of: i) adapting said input data to furtherincluding a value of a predetermined r(x_(i)) correspondingly pairedwith said predetermined x_(i); and j) constructing, via said machineprocessing unit, a value of said rational function r(x) by generating aplurality of machine processing unit signals to determine:${r(x)} = {{{r\left( x_{i} \right)} \cdot \left( {1 - \frac{{q(x)} - {q\left( x_{i} \right)}}{q(x)}} \right)} + {\frac{{p(x)} - {p\left( x_{i} \right)}}{q(x)}.}}$


15. The machine-implementable method of claim 14 wherein said rationalfunction r(x) is an approximation to a Bessel function.
 16. Themachine-implementable method of claim 14 wherein said rational functionr(x) is an approximation to an error function (ERF).
 17. Themachine-implementable method of claim 14 wherein said rational functionr(x) is an approximation to a complementary error function (ERFC). 18.The machine-implementable method of claim 14 wherein said rationalfunction r(x) is an approximation to a log gamma function (LGAMMA). 19.The machine-implementable method of claim 11 or 14 wherein all valuesare machine representations of some numerical value, said machineprocessing unit is a computer processing unit, each machinerepresentation is a binary representation operable with said computerprocessing unit, and machine-readable medium is a computer-readablemedium.
 20. A computer-readable program product having computerexecutable instructions for instructing a computer, said instructionstangibly embodying the machine-implementable method of claim
 19. 21. Acomputer-readable mathematical software routine library includingcomputer executable instructions for instructing a computer, saidinstructions tangibly embodying the is machine-implementable method ofclaim
 19. 22. The computer-readable mathematical software routinelibrary of claim 21 wherein said library is operatively associated witha software programming language.
 23. A machine for computing a propertyof a mathematically modelled physical system, the steps comprising: a)reading, via a machine processing unit, input data including a value foreach identified ordered coefficient of a first polynomial p(x)representing said property, said polynomial p(x) being expressed asp(x)=Σ(P _(j) ·x ^(j)) where j=0 to n,  a value of a quantity x, a valueof a predetermined x_(i), and a value of a predetermined p(x_(i))correspondingly paired with said predetermined x_(i); b) building, viasaid machine processing unit, a value of a second polynomial c(x) havingordered coefficients, said second polynomial c(x) being expressible as:c(x)=Σ(C _(k) ·x ^(k)) where k=0 to (n−1)  so that said first polynomialp(x) is expressible as: p(x)=p(x _(i))+{x−x _(i) }·c(x),  comprising thesteps of: i) determining, via said machine processing unit, a value foreach ordered coefficient of said second polynomial c(x) by generating aplurality of machine processing unit signals to determine each saidordered coefficient of said second polynomial c(x) from: C _(k)=Σ(P_((k+1+j)) ·x _(i) ^(j)) where j=0 to (n−1−k); ii) determining, via saidmachine processing unit, a value of said second polynomial c(x) bygenerating a plurality of machine processing unit signals to determine:c(x)=Σ(C _(k) ·x ^(k)) where k=0 to (n−1); c) constructing, via saidmachine processing unit, a value of said first polynomial p(x) bygenerating a plurality of machine processing unit signals to determine:p(x)=p(x _(i))+{x−x _(i) }·c(x);  and d) outputting, via saidmachine-processing unit, said value of the first polynomial p(x)representing said property of the mathematically modelled physicalsystem.
 24. The machine of claim 23 wherein a difference between x andx_(i) is sufficiently small to achieve a desired accuracy of a finalcomputation of said machine representation of a numerical value of saidfirst polynomial p(x).
 25. The machine of claim 24 wherein said meansfor reading said input data comprises means for reading, via saidmachine processing unit, said input data from a machine-readable medium.26. The machine of claim 25 wherein said ordered coefficients of saidsecond polynomial c(x) are computed from a mathematical expressionderivable from: C _(k)=Σ(P _((k+1+j)) ·x _(i) ^(j)) where j=0 to(n−1−k).
 27. The machine of claim 26 wherein said mathematicalexpression is a mathematical recurrence expression.
 28. The machine ofclaim 27 wherein said mathematical recurrence expression is a forwardmathematical recurrence expression.
 29. The machine of claim 27 whereinsaid mathematical recurrence expression is a backward mathematicalrecurrence expression.
 30. The machine of claim 29 further adapted toimplement said backward mathematical recurrence expression by furthercomprising: e) means for equating, via said machine processing unit, avalue of a highest ordered coefficient of said second polynomial c(x) toa value of an identified highest ordered coefficient of said firstpolynomial p(x) by generating a plurality of machine processing unitsignals to determine: C _(n−1) =P _(n);  and f) means for computing, viasaid machine processing unit, a value for each lower ordered coefficientof said second polynomial c(x) by generating a plurality of machineprocessing unit signals to determine: C _(k+1) =x _(i) ·C _(k) +P _(k)where k=(n−1) to
 1. 31. The machine of claim 30 wherein saidpredetermined x_(i) is selected from a set of predetermined values ofx_(i).
 32. The machine of claim 30 wherein said predetermined x_(i) is aclosest member of said set to said identified x.
 33. The machine ofclaim 32 wherein the determining means for determining a value of saidsecond polynomial c(x) is computed by using Homer's Rule.
 34. Themachine of claim 33 for determining a value of a denominator polynomialq(x) having identified ordered denominator coefficients, saiddenominator polynomial q(x) being expressible as: q(x)=Σ(Q _(j) ·x _(j))where j=0 to m, comprising further steps of: g) adapting said input datato further include a value for each identified ordered denominatorcoefficient of said denominator polynomial q(x), and a value of apredetermined q(x_(i)) correspondingly paired with said predeterminedx_(i); and h) determining, via said machine processing unit, a value ofanother polynomial d(x) having ordered denominator coefficients, saidanother polynomial d(x) being expressible as: d(x)=Σ(D _(k) ·x ^(k))where k=0 to (m−1)  so that said denominator polynomial q(x) isexpressible as: q(x)=q(x _(i))+{x−x _(i) }·d(x),  and a value for thesaid denominator polynomial is resolved.
 35. The machine of claim 34wherein the first polynomial p(x) is a numerator polynomial p(x), andp(x)−p(x_(i)) is computed, and p(x_(i)) is not read.
 36. The machine ofclaim 35 for determining a value of a rational function r(x) beingexpressible as a quotient of said numerator polynomial p(x) and saiddenominator polynomial q(x) expressed as${{r(x)} = \frac{p(x)}{q(x)}},$

comprising further steps of: i) adapting said input data to furtherincluding a value of a predetermined r(x_(i)) correspondingly pairedwith said predetermined x_(i); and j) constructing, via said machineprocessing unit, a value of said rational function r(x) by generating aplurality of machine processing unit signals to determine:${r(x)} = {{{r\left( x_{i} \right)} \cdot \left( {1 - \frac{{q(x)} - {q\left( x_{i} \right)}}{q(x)}} \right)} + {\frac{{p(x)} - {p\left( x_{i} \right)}}{p(x)}.}}$


37. The machine of claim 36 wherein said rational function is anapproximation to a Bessel function.
 38. The machine of claim 36 whereinsaid rational function is an approximation to an error function (ERF).39. The machine of claim 36 wherein said rational function is anapproximation to a complementary error function (ERFC).
 40. The machineof claim 36 wherein said rational function is an approximation to a loggamma function (LGAMMA).
 41. The machine of claim 33 or 36 wherein allvalues are machine representations of some numerical value, said machineprocessing unit is a computer processing unit, each machinerepresentation is a binary representation operable with said computerprocessing unit, and said machine-readable medium is a computer-readablemedium.
 42. A machine having a computer-readable program product havingcomputer executable instructions for instructing a computer to embodythe machine of claim
 41. 43. A machine having a computer-readablemathematical software routine library including computer executableinstructions for instructing a computer to embody the machine of claim41.
 44. A machine having the computer-readable mathematical softwareroutine library of claim 43 wherein said library is operativelyassociated with a software programming language.